DESY 09-159 
SFB/CPP-09-90 



Continuum-limit scaling of overlap fermions 
as valence quarks 

KRZYSZTOF CICHyE 

Adam Mickiewicz University, Faculty of Physics, 
^ I Umultowska 85, 6I-6I4 Poznan, Poland 

^ ; GREGORIO HERDOIZA, KARL JANSEN 

NIC, DESY, Platanenallee 6, D-15738 Zeuthen, Germany 

o 



in 



Oh: 



o 



X 



Abstract 



We present the results of a mixed action approach, employing dy- 
namical twisted mass fermions in the sea sector and overlap valence 
\ fermions, with the aim of testing the continuum limit scaling behaviour 

of physical quantities, taking the pion decay constant as an example. 
To render the computations practical, we impose for this purpose a 
\ fixed finite volume with lattice size L ~ 1.3 fm. We also briefly review 

the techniques we have used to deal with overlap fermions. 



> ■ PACS numbers: ll.15.Ha, 12.38.Gc 

O 1 Introduction 

d 



Lattice QCD is considered to be the most effective way of studying the 
non-perturbative aspects of the field theory of strong interactions, Quan- 
tum Chromodynamics (QCD). It is a regularization of QCD, which consists 
in discretizing the relevant degrees of freedom by putting them on a four 
dimensional lattice with lattice spacing a, its inverse being the ultraviolet 
5^ \ cutoff of the theory. With the present generation of supercomputers, large 

scale dynamical simulations are performed by a number of collaborations us- 
ing different types of gauge and fermionic lattice actions. However, for some 
classes of actions, such as overlap fermions, which exactly preserve chiral 
symmetry, such simulations are extremely demanding. A promising alter- 
native for fully dynamical simulations with overlap fermions is the mixed 



*Presented at the 49. Cracow School of Theoretical Physics, Zakopane, Poland, 31 
May - 10 June 2009. 

^For instance, the JLQCD and TWQCD collaborations are simulating dynamical over- 
lap fermions, but in their case the global topological sector is kept fixed [T]. 
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action approach, which consists in using a computationally faster formula- 
tion for the fermions in the sea, such as Wilson twisted mass fermions, while 
using overlap fermions in the valence sector. In this way, one avoids the most 
computer-intensive part of a dynamical simulation with overlap fermions, but 
at the same time one profits from having such fermions in the valence sector. 
In the case of overlap fermions, its exact chiral symmetry gives the benefit 
of simplifying the operator mixing problem, something which is essential in 
several types of lattice computations, like the determination of the kaon bag 
parameter Bk- 

In this paper, we present the results of a mixed action approach with 
overlap fermions in the valence sector and Wilson twisted mass fermions in 
the sea sectoiH. In Section 2 we briefly review the overlap formulation and the 
techniques necessary to effectively apply it. Section 3 gives the description 
of our setup and in Section 4 we present the continuum-limit scahng test of 
overlap fermions and we discuss the results. Section 5 concludes. 

2 A brief review of overlap fermions 

Since this contribution is written for a more general audience, we provide 
here a very short review of overlap fermions. 

2.1 The need for overlap fermions 

A very general problem of lattice field theory with fermions is the doubling 
problem. With a naive fermion discretization, instead of one fermion in the 
continuum limit, one has 16 fermions in 4 dimensions of spacetime. As was 
originally proposed by Wilson [1], the doubling problem can be solved by 
using the following Wilson-Dirac operator: 

Dw{m) = ^ (7^(V; + V^) - arV;V^) + m, (1) 

where V* (Vfj,) is the backward (forward) lattice derivative, r - Wilson pa- 
rameter, m - bare quark mass. However, the Wilson term in the above Dirac 
operator explicitly breaks chiral symmetry, which is in accordance with a 
general theorem proved by Nielsen and Ninomiya [S] back in 1981 that it is 
impossible to have at the same time for a Dirac operator D: locality, trans- 
lational invariance, no doublers and chiral symmetry, i.e. {D(p),75} = 0. 
Since Wilson's proposal, much of the effort went into finding a lattice theory 
without doublers which preserves the largest possible number of symmetries, 

^For an account of earlier stages of the project, see [2], [S]. 
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and at the same time reaches the continuum hmit as fast as possible (which 
practically means the absence of 0{a) leading cutoff dependence). 

In 1982 (i.e. only one year after establishing the Nielsen- Ninomiya the- 
orem), it was shown by Ginsparg and Wilson [6] that a remnant of chiral 
symmetry is present on the lattice without the doubler modes, if the cor- 
responding Dirac operator D obeys an equation now called the Ginsparg- 
Wilson relation: 

-f,D + D75 = aD^.D. (2) 

However, the solutions to this equation were not known for many years, until 
a particularly simple form of a Dirac operator obeying the Ginsparg- Wilson 
relation was given by Neuberger ^ in 1997. Neuberger's discretization is 
now usually referred to as overlap fermion^ and the (massless) overlap Dirac 
operator is of the form: 

Do,(0) = -(1-A(AU)-V2), (3) 

where: 

A = l-aDwiO). (4) 

It was also found that the overlap Dirac operator is local under very general 
conditions [H]. Moreover, in 1998 Liischer [TU] found that the Ginsparg- 
Wilson relation leads to a non-standard realization of lattice chiral symmetry. 
The action is invariant under: 

^^e''^'-"^)^, (5) 

^^^e'e75(i-^), (6) 

Although it is not the standard form of chiral symmetry, it still correctly 
reproduces the anomaly and protects the fermions from additive mass renor- 
malization and 0{a) lattice artifacts. The non-standard realization of chiral 
symmetry means also that the conditions of the Nielsen- Ninomiya theorem 
do not apply and one can have chiral symmetry (which becomes standard 
chiral symmetry in the continuum limit) without the doublers. 

In order to simulate massive overlap quarks, one uses the following form 
of the Dirac operator: 

D^m) = {1 - DUO) + m, (7) 

where m is the bare overlap quark mass and Dov{0) the massless overlap 
Dirac operator of eq. ([3]). 

■^For a review of overlap fermions see e.g. [5]. 
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2.2 Techniques to effectively deal with overlap fermions 



The main disadvantage of overlap fermions is that they are much more costly 
to simulate - by a factor of 30-120 in comparison with maximally twisted 
mass fermions [TT]. Moreover, this factor increases when approaching the 
physical pion mass. Therefore, it was essential to develop techniques to 
effectively deal with overlap fermions. The aim of this subsection is to briefly 
review them. 

Computation of the overlap operator. The first important thing is to 
effectively calculate the overlap Dirac operator itself. This is non-trivial, since 
it involves the inverse of the square root of a matrix. There are several ways 
to do this, including polynomial approximations, Lanczos based methods 
and partial fraction expansion. For overviews of these methods, see e.g. 
|12j . |13j . The method that we have chosen is a Chebyshev polynomial 
approximation, which consists in constructing a polynomial P„^e(x) of degree 
n, which has an exponential convergence rate in the interval x G [e, 1]. The 
advantages of using this approximation are the well-controlled exponential 
fit accuracy and the possibility of having numerically very stable recursion 
relations which allows for large degrees of the polynomial. To ensure that 
the Ginsparg- Wilson relation (for massless Dirac operator) is fulfilled to a 
very high precision, the degree of Chebyshev polynomial n has to satisfy the 
following condition: 

||X - P^,M^A)A^AP^,,{A^A)X\\y\\X\\^ < (8) 

where ^ has to be a very small number, typically set to 10~^^ to achieve a 
compromise between good quality of approximation and its cost. The degree 
of the polynomial depends on the condition number of the matrix A^A, i.e. 
the ratio of the highest to the lowest eigenvalue. The lowest eigenvalue can 
be a very small number and hence the condition number can be prohibitively 
large, if one constructs the approximation on the interval [e, 1], with e being 
the lowest eigenvalue. Fortunately, one can do much better with the following 
method. 

Eigenvalue deflation. To achieve a considerably smaller degree of Cheby- 
shev polynomial, one should calculate a certain numbei0 of the lowest eigen- 
values and eigenvectors of A'^A and project them out of this matrix. In this 
way, the Chebyshev approximation is constructed on the interval [e, 1] with e 
equal to the highest of the computed eigenvalues. To illustrate how low the 
lowest eigenvalues of the matrix A^A can go, we plotted in Fig. [1] the five 
lowest eigenvalues and the highest eigenvalue for four gauge field ensembles, 

^The exact number has to be found empirically and tends to increase with volume. 
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Figure 1: 5 lowest eigenvalues and the highest eigenvalue for various gauge 
field ensembles. The lattice spacing is a ~ 0.084 fm for upper plots, a ~ 0.066 
fm for bottom left and a ^ 0.053 fm for bottom right plot. 



each with 200 configurations. The general dependence that can be deduced 
from these plots is that increasing lattice spacing (decreasing f3) moves the 
spectrum down (eigenvalues in lattice units tend to become smaller) and 
increases the probability of having very low eigenvalues. 
HYP smearing of gauge fields. Another way to lower the condition num- 
ber of the matrix A'^A and thus the degree of the Chebyshev polynomial is to 
perform one iteration of HYP smearing on the gauge fields. This technique 
was introduced by A. Hasenfratz and F. Knechtli [14J and allows to achieve 
much better convergence of the solver due to improved chiral propertied. In 
comparison with other link fattening methods (e.g. APE smearing), HYP 
smearing is beheved to preserve better the short-distance quantities, because 
it mixes links from hypercubes attached only to the original link. 
SUMR solver with adaptive precision and multiple mass capability. 
Having constructed the Dirac operator (with Chebyshev approximation), to 

^Hasenfratz and Knechtli [H] remark that fat hnks lead to an order of magnitude 
improvement in convergence. 
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find the propagator one has to solve the equation: 

Dov{m)tjj = 1], (9) 

where ijj is the propagator and rj is the source - a vector whose choice will be 
commented on below. To effectively solve this equation, one has to choose 
the most appropriate solver. Chiarappa et al. [H] found that in the case 
of (quenched) overlap and small volume (12^ and 16^), the most effective 
solver is the chiral conjugate gradient algorithm, with the SUMR solver just 
behincl^. However, the former algorithm can only be used for exact overlap 
operator, which means that the polynomial approximation would lead to 
some corrections that would have to be explicitly calculated. This makes the 
latter algorithm more attractive and we decided to use it. To improve its 
performance, we have also used adaptive precision and multiple masses. The 
former means that the degree of the Chebyshev polynomial is adapted to 
what accuracy is actually needed in the present iteration step. In practice, 
this means that when the solver is moving towards the requested precision, 
the accuracy of approximation can be substantially decreasedll, which saves a 
factor of around 2 in inversion time. Multiple mass capability of an algorithm 
means that for the cost of one inversion for the smallest bare quark mass, 
one can obtain the solution also for heavier quark masses for practically no 
additional cost. Since the dependence on the quark mass is central in the 
present project, the use of a multiple mass variant of the SUMR solver was 
absolutely essential. 

Stochastic sources. Another important aspect of solving eq. iQ is the 
choice of the source rj. The most obvious choice is the point source, which 
means that the vector rj is chosen to be 1 at one space-time point and spin- 
color component and otherwise. This leads to 12 inversions for each gauge 
configuration, one per spin-color component. However, one can do much 
better when it comes to statistical error on the pion mass and especially the 
pion decay constant, with the use of stochastic sources. To keep the signal-to- 
noise ratio high enough, one has to use timeslice sources, i.e. sources that are 
non-zero for all spatial points on a given time-slice and for a given spin. This 
requires 4 inversions per gauge configuration, one per spin component. More- 
over, the noise can be further reduced by using the one-end trick, introduced 
in [T3]. A reasonable choice of stochastic sources for mesonic correlators 
are the Z{2) sources - the random numbers are of the form (±1 ± t)/V2. 

•^SUMR = Shifted Unitary Minimal Residual. 

^For example, if the degree of Chebyshev polynomial at the start of inversion is typically 
(for our parameters) of order 150-200, in the final iterations it can go down typically to 
30-40 with adaptive precision. 



6 



Empirical observations show that this method reduces the statistical error 
on the pion decay constant by a factor of approximately 2 for the smallest 
quark masses that we consider, which means that the number of inversions 
to achieve the same statistical error can be even four times smaller than with 
point sources. 

3 Setup 

To perform the scaling test of overlap fermions we fix the volume to L ^ 1.3 
fm and use the following ensembles of dynamical Nf = 2 maximally twisted 
mas^ [I7j, [IB] configurations, generated by the European Twisted Mass 
Collaboration (ETMC) [IH], [20]: 

• P = 3.9, = 16^ X 32, a ^ 0.084 fm, a/i = 0.004, k = 0.160856, 

• P = 4.05, V = 20^ X 40, a ^ 0.066 fm, afx = 0.003, k = 0.157010, 

• P = 4.2, V = 24^ X 48, a ^ 0.053 fm, a/i = 0.002, k = 0.154073. 

The valence quarks are overlap fermions and the overlap quark mass was 
chosen to vary from the unitary light quark mass up to the physical strange 
quark mass. In total we have 20 quark masses, which allows for a precise 
determination of the quark mass dependence of the pion mass and decay 
constant. 

In addition, we also perform a tree- level scaling test, for which we fix 
Nm = 0.5 (which is the equivalent of fixing volume), where A'^ is the number 
of lattice points in spatial directions and go from = 4 to = 64. Thus, 
the change in introduces the scaling towards the continuum limit, which 
corresponds to A^ — > cxd. The temporal extent was set to be 64 times larger 
than the spatial extent for each value of A^, which makes it possible to extract 
the relevant quantities without any contaminations from the excited states. 
For the details of this test and expressions for the tree-level quark propagators 
and correlation functions, see ^21j, ^22j. 

^Twisted mass (TM) Dirac operator is defined by: Dtm = £'vK(m)l/ + 1/^75 r^, where: 
m - untwisted quark mass, jj, - twisted quark mass, If and act in flavour space. For a 
review of twisted mass fermions, see |16) . 
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Figure 2: The relative cutoff effects of the pseudoscalar correlator at a fixed 
physical distance, pseudoscalar mass and decay constant. 

4 Results 

4.1 Tree- level test 

We investigated the relative (with reference to the continuum-limit value) 
cut-off effects of three observable^: pseudoscalar (PS) correlation function 
N^Cps at a fixed physical distance t/N = 4, pseudoscalar meson mass Nmps 
and pseudoscalar decay constant N fps- The results are presented in Fig. [21 
As expected, all observables show only 0{a^) scaling violations, since over- 
lap fermions are (9(a)-improved. We also observe that the smallest lattice 
artifacts are observed for the pseudoscalar mass and the largest for the cor- 
relation function itself. 

4.2 The interacting case — matching the pion mass 

In the interacting case, we are interested in the continuum limit scaling of 
the pion decay constant for three reference values of the pion mas^ - the sea 
quark pion mass, an intermediate mass romj^ ~ 1 and a heavy mass (around 

^AU of the observables are multiplied by an appropriate power of N to make them 
dimensionless. 

^"The pion mass and decay constant were obtained from the pseudoscalar correlation 
function. Thus, pion means here the light pseudoscalar meson. 
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Figure 3: Matching mass for ensemble 16^ x 32, a ~ 0.084 fm {(3 = 3.9), 
a/i = 0.004, K = 0.160856. The matching procedure gives am = 0.007(1). 



the strange quark region) rQin.,^ ^ 1.50 

The matching quark mass m is defined, for each ensemble, by the condi- 
tion: 

mr(m)=m™(/i), (10) 

where and m™ are the mixed action and unitary pion masses, respec- 
tively. 

Fig. [3] shows the procedure of matching for the ensemble at /? = 3.9. 
The "overlap" curve shows the dependence of the pion mass on bare overlap 
quark mass. As the plot indicates, when the overlap quark mass equals 0.007 
(with an uncertainty of approx. 0.001), the TM pion and overlap pion have 
equal masses, corresponding in infinite volume to ca. 300 MeV. An analogous 
procedure was applied for other ensembles. 

One should add here a word of caution. For the situation of rather small 
lattice extents of L ^ 1.3 fm, the chiral zero- modes of the overlap operator 
can play an important role. Note that since we use here a mixed action 
approach, these zero-modes are not compensated for by the fermion deter- 
minant. The infiuence of the interplay between the finite box size, the quark 
mass and the mixed action on the effects of the zero-modes for physical ob- 
servables will be discussed elsewhere. 

These two higher reference masses correspond to the partially quenched setup. 
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Figure 4: Continuum limit scaling of the (overlap) pion decay constant for 
fixed volume L ^ 1.3 fm. 



4.3 Continuum limit scaling of the pion decay constant 

Fig. m shows the continuum limit scaling of to/tt for three choices of the 
pion mass. The lower curve corresponds to the sea quark pion mass, i.e. 
the bare overlap quark mass is set to the relevant m for each ensemble. For 
the intermediate and upper curve we fix the pion mass rom^ to 1.0 and 1.5, 
respectively, i.e. we take such overlap quark mass that leads to the chosen 
value of the pion mass for each ensemble. 

We expect O(a^) scaling violations and hence we plot the decay constant 
against the lattice spacing squared. 

Indeed, for all three cases, we observe good scaling behaviour towards 
the continuum limit. The fact that the unitarity violations, proper to a 
mixed action approach, do not spoil the scaling of the pion decay constant 
is reassuring. 



5 Conclusion and prospects 

In this paper, we briefly reviewed the overlap discretization of the fermio- 
nic action and the techniques we have used to deal with overlap fermions. 
We presented the continuum limit scaling test of the pion decay constant at 
tree-level and in the interacting case, using three values of the lattice spac- 
ing coming from dynamical twisted mass configurations at roughly matched 



10 



physical box length of L ^ 1.3 fm and with three reference values of the pion 
mass. 

We observe very good continuum limit scaling properties for all cases, 
i.e. at tree-level and in the interacting case for pion masses corresponding 
to light (with overlap quark mass leading to the same pion mass as the sea 
quark), intermediate and somewhat heavier quarks. 

In the next stage of the project, we plan to extend the scaling test to 
other observables (e.g. the nucleon mass) and we aim at a detailed study 
of the particular mixed action setup that we use, including fits of Partially 
Quenched Chiral Perturbation Theory formulae. We will also investigate the 
possibility of computing observables for which chiral symmetry is crucial. 
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